
module Distribution.Hurst
        ( hurstExp
        )
        where

import Numeric.GSL.Fitting.Linear
import Data.Packed.Vector

mean ::  Fractional a => [a] -> a
mean xs = (sum xs)/(fromIntegral $ length xs)

meanAdjusted ::  Num b => b -> [b] -> [b]
meanAdjusted m xs = map (\x -> x - m) xs

cumulativeDeviated ::  Num a => [a] -> [a]
cumulativeDeviated xs = scanl1 (\x y -> x+y) xs

rangeR ::  (Ord a, Num a) => [a] -> a
rangeR xs = (maximum xs) - (minimum xs)

-- mean-adjusted is assumed
stddev ::  Floating a => [a] -> a
stddev xs = sqrt (sqrsums/n)
        where   n       = fromIntegral $ length xs
                sqrsums = sum $ map (\x -> x**2) xs

rescaledRange ::  (Ord a, Floating a) => [a] -> a
rescaledRange xs = r/s
        where   r       = rangeR cumm
                m       = mean xs
                madj    = meanAdjusted m xs
                cumm    = cumulativeDeviated madj
                s       = stddev madj

vectorize :: [(Double, Double)] -> (Vector Double, Vector Double)
vectorize = buildVectors . unzip
        where  buildVectors (x, y)     = (fromList x, fromList y)

remapToRs ::  [(t, [Double])] -> [(t, Double)]
remapToRs = map (\(n, xs)->(n, rescaledRange xs))

hurstExp ::  [(Double, [Double])] -> Double
hurstExp dts = k
        where (_, k, _, _, _, _) = ((uncurry linear) . vectorize . remapToRs) dts

